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CHAPTER  ONE 
INTRODUCTION 

The  ability  to  calculate  the  fluid  dynamics  that  result  from  flame 
spreading  and  pressure  wave  formation  in  a bed  of  highly  loaded,  small 
grain,  solid  propellant  is  the  motivation  for  the  work  reported  here. 
There  have  been  previous  efforts  at  the  University  of  Illinois,  under 
the  direction  of  Dr.  Herman  Krier,  to  analyze  such  transient  flow 
problems.  These  include  the  work  reported  by  Van  Tassell  and  Krier 
[1],  Krier  and  Gokhale  [2J,  Krier,  Rajan,  and  Van  Tassell  [3],  Dimitstein 
[4],  Krier,  Dimitstein,  and  Gokhale  [5],  and  Krier,  Gokhale,  and 
Hughes  [6].  This  considerable  background  proved  quite  useful  and  indeed 
essential  in  allowing  this  detailed  work  to  be  completed  in  a relatively 
short  period  of  time. 

There  have  been  several  different  approaches  in  developing  the 
appropriate  governing  conservation  equations  to  describe  the  physics  of 
the  flow.  Although  it  is  difficult  to  categorize  all  the  models 
published,  one  of  the  main  differences  stems  from  the  assumptions  that 

(a)  the  equations  of  motion  be  derived  based  on  the  center  of  mass  of 
each  phase,  versus  (b)  those  derived  based  on  the  center  of  mass  of  the 
mixture.  In  addition  to  that  distinction,  there  are  models  that  (c) 
attempt  to  define  rigorous  average-flow  properties  in  which  one  phase 
does  not  constitute  a continuum. 

A brief  review  was  recently  prepared  by  Krier  [7]  which  discusses 
the  concepts  resulting  in  equations  based  on  assumption  (a)  versus 
those  obtained  by  assumption  (b) . Derivations  utilizing  assumption  (a) 
have  been  termed  the  "separated-flow"  approach,  those  using  assumption 

(b)  are  called  the  "continuum-mixture"  approach,  and  tliose  using 
assumption  (c)  might  come  under  the  broad  definition  of  an  "averaging" 
approach. 

The  work  reported  here  is  an  attempt;  to  first,  derive  the  con- 
servation equations  using  the  concept  of  separated  flow  (as  defined  in 
the  text  by  Wallis  [8])  and  second,  to  develop  an  accurate  solution 
technique  to  solve  the  problem  of  unsteady  flows  through  packed  beds  of 
propellant.  The  derivation  of  the  conservation  equations  is  presented 
in  Appendix  A.  The  equations  derived  here  are  compared  with  those  of 
other  investigators  and  discussed  in  Chapter  II. 


1 


2 

This  work  also  includes  some  modeling  that  will  allow  the  calcu- 
lations to  be  utilized  at  high  gas  pressures.  As  discussed  in  the 
text  that  follows,  such  items  as  (a)  an  accurate  non-ideal  high 
pressure  equation  of  state,  (b)  a better  value  for  the  mixture 
sound  speed,  (c)  a model  for  the  axial  particle  stress,  and  (d)  a 
temperature  dependent  viscosity  were  included.  With  all  of  the 
above  inputs  it  will  be  shown  that  this  code  provides  reasonable  and 
consistant  pressure  wave  and  flame  front  histories  over  widely  varied 
flow  regimes. 

Def lagration-to-Detonation  Transition 

Studies  of  the  spontaneous  def lagration-to-detonation  transition 

(DDT)  in  gases  indicate  that  in  an  accelerated  deflagration  a shock 

front  runs  ahead  of  the  flame.  The  present  work  is  concerned  with 

such  dynamic  behavior  in  granular  propellant,  with  an  attempt  to 

elucidate  the  details  of  propagation  of  the  flame  (ignition)  front, 

pressure  fronts  and  possible  shock  fronts  before  development  of  steady- 

35 

state  detonation.  The  experimental  work  of  Bernecker  and  Price  with 
granular  explosives  and  the  work  of  Gipson  and  Macek^^  with  DDT  in 
solid  explosives  indicate  clearly  that  in  the  study  of  buildup  to 
detonation  it  is  important  to  know  the  relationship  of  the  pressure 
(shock  or  compression)  fronts  to  the  flame  front. 

In  the  problem  of  concern  here  a fraction  of  the  propellant  grains 
are  assumed  ignited  at  one  end  of  a closed  chamber.  The  mass  generated 
in  the  ignition  region  accelerates  the  flame  forward,  with  the  region 
behind  the  flame  rapidly  increasing  in  pressure  as  more  and  more  mass 
is  being  generated  due  to  the  pressure-sensitive  burning  rate.  Cal- 
culations then  show  that  a steep  pressure  front  begins  to  accelerate 
and  attempts  to  overtake  the  flame  front.  When  this  happens  there  is 
an  almost  abrupt  increase  in  the  flame  speed.  A detonation  can  be  said 
to  have  started  when  the  pressure  front  proceeds  the  flame  front  and 
they  both  attain  about  the  same  high  speed.  We  are  of  course  limited 
by  our  numerical  differencing  scheme  to  indicate  a true  shock  front 
proceeding  the  flame.  But  with  a rapid  transition  of  the  flame  front 
in  conjunction  with  the  interaction  of  the  pressure  fronts,  one  can 
suppose  that  DDT  can  occur.  A discussion  of  a criterion  for  DDT  is 
presented  later. 
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CHAPTER  TWO 

THE  CONSERVATION  EQUATIONS  FOR  ONE -DIMENSIONAL  I 

UNSTEADY  TWO-PHASE  FLOW 

2.1  Introduction 

' I 

The  one-dimensional  conservation  equations  for  a two-phase  re- 
active mixture  have  been  derived  by  most  investigators  through  the 
concept  of  separated  flow  as  defined  in  the  text  by  Wallis  [8  ]. 

This  approach  considers  two  distinct  flows,  each  through  a separate  j 

control  volume,  such  that  the  sum  of  the  volumes  represents  an  average  i 

mixture  volume  and  the  sums  of  the  properties  of  each  flow  represent  ; 

average  mixture  properties  of  the  fluid.  With  this  approach,  separate 
equations  for  continuity,  momentum,  and  energy  are  written  for  each 
phase  and  are  solved  with  a state  equation  to  describe  the  overall 
flow.  Panton  [10]  has  provided  the  rigorous  formulation  with  the 
important  assumptions  required  to  express  these  field  balance  equations. 

There  are,  of  course,  many  other  analyses  that  have  arrived  at  similar  J 

forms  of  the  equations.  Some  of  these  would  include  the  work  by  Ander- 
son and  Jackson  [ll],  Pigford  and  Baron  [12],  Green  and  Naghdi  [13], 

Crowe  [14],  Solbrig,  Mortenson,  Lyczkowski  [15],  Hughes  [l6],  Kuo  and 
Summerfield  [17] , Gough  [18],  Nigmatulin  [19],  Gidaspow  [20],  and 
Zuber  [21].  Not  all  of  these  investigations  have  treated  the  two  energy 
equations,  as  is  done  here, 

2.2  Assumptions 

The  important  assumptions  upon  which  this  separated  flow  model  is 
based  are  listed  below. 

1)  The  two  phases  are  assumed  interdispersed , but  separate, 
coupled  only  by  the  appropriate  interaction  terms. 

2)  Each  phase  is  a continuum  and  therefore  derivatives  are  unique- 
ly defined. 

3)  The  total  cross  sectional  area  is  the  sum  of  the  gas  and  solid 
flow  areas.  (The  problem  is  quasi-one-dimensional) . 

4)  Solid  particles  are  large  compared  to  molecules  and  hence 
do  not  contribute  to  the  total  pressure  of  the  mixture. 

5)  When  combustion  of  particles  causes  mass  transfer  between  the 

phases,  the  solid  phase  always  loses  mass  while  the  gas  phase 
gains  it,  i.e.  , ^ 0, 
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6)  The  diameter  of  a typical  particle  can  be  thought  of  as 
the  average  diameter  of  the  total  number  of  particles. 

7)  All  gases  obey  the  non-ideal  Noble-Abel  equation  of  state 
with  a variable  co-volume. 

8)  Heat  transfer  to  the  particles  due  to  conduction  and  radiation 
is  neglected. 

9)  The  density  of  the  solid  phase  is  constant. 

10)  The  equations  are  laminar  in  the  sense  that  turbulent  flow 

due  to  the  two-phase  nature  of  the  fluid  has  been  averaged  out. 
(See  Panton  [10]  for  additional  discussions  regarding  this 
assumption. ) 

11)  The  gases  are  inviscid  except  for  their  action  on  the  particles 
through  the  drag  term. 

12)  The  fluid  properties  and  are  assumed  constant,  but  gas 
viscosity  and  conductivity  are  generally  known  functions  of 
temperature. 

2.3  The  Conservation  Field  Balance  Equations 

Some  of  the  details  in  the  derivation  of  the  six  field  balance 
equations  have  been  presented  in  Appendix  A.  In  order  to  specifically 
describe  the  state  of  a two-phase  fluid  one  needs  equations  to  solve 


for  the  following  variables: 


P . u , 

g g 


g 


P , u , and  T . 
g P P 


The 


state  equation  is  added  to  the  field  balance  equations  to  obtain  the 
seven  required  equations.  Therefore  two  continuity  equations,  two 
momentum  equations,  and  two  energy  equations  will  be  necessary.  These 
equations  are  listed  below  in  operator  form.* 

Gas  Continuity 


3p 

W 


1 


3u 

^1  ^ ■ “g  air 


^ r 

S 3x 


g 


Solid  Continuity 


3P2  _ ^ ^2^  3S  . 

~ " ^2  3x  '^p  3x  S 3x  g 


3t 

Gas  Momentum 
3u 

— g.  = - 

3t 


A ->n  (U  -U  ) 

_B.  . i H 


g 


Pj  3x 


'1 


r - — 

g Pi 


(2.1) 


(2.2) 


(2.3) 


See  list  of  nomenclature  for  definitions  of  all 
symbols  in  equations  (2.1)  - (2.6). 


5 


Solid  Momentum 


^ = - 
3t 


,,  3t  tt 

E.  - E + L 


p 3x 


Gas  Energy 
3E 

= . 
3t 


3E 

1 — S- 

3x 


r 

Pi 


(E. 


3x 


u P 


PlS 


3S 

3x 


Pi 


3(^ 

3x 


Pi 


n ^ 


r u 

-S.  r -S.  _ 


3u 

_& 

3x 


(2.4) 


(2.5 


2 _ 


u u 
g P 


u 

+ -^)- 


-(u  -u  )-- 
g p Pi 


Solid  Energy 

3E  3E  .vU  3t  r u * 

3t  p 3x  Pi  9x  P,  ^ P p p_ 


(2.6) 


Here,  specifically 


Pi  = <^Pg 


A 


r = 

g 


c n ^ 1 2 

E = E + ::7  U 

go  g 2 g 

C C ^ 1 2 

E = E + •=•  u 

Po  P 2 p 


P2  = (l-fti)  P 


1 - ({)  = ^ 

E = C T 
g V g 


E = C T 
P P g 


The  subscript  "g"  denotes  gas  and  "p"  denotes  particle.  The  subscript 
"o"  stands  for  (totcl)  stagnation  conditions. 


The  seventh  equatiw  is  the  Noble-Abel  state  equation 

p R T 

p = g g g 
g 1 - p B 

6 V 


(2.7) 


where  is  the  covolume  which  is  a function  of  gas  density,  so  that 
at  extreme  pressures,  an  appropriate  non-ideal  equation  is  satisfied. 

There  is  currently  concern  whether  the  quasi-linear,  one- 
dimensional, separated  flow  equations  listed  above,  or  those  derived 
from  other  formulations,  constitute  an  hyperbolic  system.  A recent 
examination  of  such  systems  with  regard  to  this  question  is  carried 
out  by  Hughes  [6].  It  is  the  contention  here  that  the  system  of 
equations  must  be  hyperbolic  to  provide  stable  and  meaningful  results. 
It  is  also  felt  that  knowledge  of  the  nature  of  the  roots  (character- 
istics) of  the  equations  can  provide  much  information  concerning  the 
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boundary  conditions,  initial  conditions  and  other  auxiliary  information 
required  to  solve  the  system.  Therefore  changes  in  the  form  of  the 
six  field  balance  conservation  equations  (in  particular  the  two  momen- 
tum balances)  may  be  required  in  order  to  provide  a well-posed,  hyper- 
bolic system. 

Another  reason  for  alternate  forms  of  these  balance  equations 
arises  from  the  derivation  of  the  solid  phase  energy  equation.  There 
can  be  some  ambiguity  concerning  the  definition  of  the  appropriate 
total  enthalpy  of  the  particles.  Equation  2.6  above  assumes 

Hp^  = Ep  . Y ^2.8) 

instead  of  the  usual  definition 


Since  the  meaning  of  a particle  pressure,  P^,  is  unclear,  some  investi- 
gators, like  Rudinger  [27]  have  assumed  that 


P u 

H = E + -£-  + -4 

Po  P "p  2 


(2.10) 


is  more  appropriate.  All  of  these  definitions  for  H imply  a possible 

Po 

interchange  between  the  particle  internal  energy  and  kinetic  energy. 
Another  possible  assumption  that  does  not  allow  for  such  an  energy  trans- 
fer is 


H = E 
P P 


(2.11) 


only.  Assuming  that  the  total  stagnation  enthalpy  is  only  the  internal 
energy  appears  plausible  because  first,  the  solid  phase  is  incompressible 
and  second,  it  does  not  constitute  a gas.  Therefore  a change  in  kinetic 
energy  should  not  necessarily  imply  a change  in  temperature. 

Particle  energy  equations  based  on  the  three  assumptions  listed 
above  (equations  2.8,  2.10,  and  2.11)  are  derived  in  Appendix  A.  Final 
forms  for  the  latter  two  assumptions  are 
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(l-<}i)P  3u  u P - (l-d))u  3P 

g P _ P g p g 


(l-(())u  3t  r u P 4> 

E ^ + _g.  [E  + + -E 

P2  P2  P P SPp  9x  P2 


for  the  assumption  in  equation  2.10,  and  the  following. 


uj-  VI-  , u ru 

^ . U 5-E  . J [E  . e‘=''™  . 4-]  . -E  . 


(2.12) 


(2.13) 


for  Hp  = Ep.  The  systems  based  on  all  three  variations  were  studied 
and  di?ferences  are  reported  in  Chapter  4. 

It  has  been  duly  observed  that  several  investigators,  in  particular 
Kuo  ej^  al . , [17,22]  and  Gough  [18],  do  not  write  a solid  phase  energy 
balance  s>'mmetric  with  the  gas  phase  energy  balance,  like  equation  2.6 
above.  Instead,  these  investigators  have  decided  to  simply  solve  the 
unsteady  heat  conduction  equation  for  the  solid  particle  and  monitor 
the  heat  gained  from  the  gas  in  the  form  of  convective  heat  transfer. 
This  approach  was  not  deemed  acceptable  for  the  work  done  here.  It  is 
considered  fundamentally  necessary  that  when  adding  the  energy  equations 
2.5  and  2.6,  the  classical  mixture  energy  equation  is  obtained  as 
required  in  the  discussion  by  Wallis  [8].  This  leads  to  the  utilization 
of  an  overall  or  "lumped"  parameter  for  particle  energy  (temperature) 
which  requires  an  appropriate  ignition  criteria.  These  discrepancies 
between  models  are  summarized  in  Table  2-2. 

2.4  Comparison  with  Other  Models 


Table  2.1  lists  both  the  gas  and  particle  phase  momentum  equations 
as  reported  by  six  groups  of  investigators.  In  preparing  the  table 
some  work  was  required  because  several  of  the  equations  were  reported 
only  in  conservative  form.  As  demonstrated  in  Appendix  A,  the  con- 
tinuity equation  must  be  utilized  to  transform  the  momentum  field 
balance  equations  into  "operator"  form. 

The  left  hand  side  (L.H.S.)  of  both  field  balance  equations  are 
identical  for  all  models,  as  expected.  The  two  momentum  equations  used 
here  generally  agree  well  with  those  of  other  investigators  using  the 
separated  flow-concept.  There  are  some  differences,  however. 
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In  the  form  used  by  Krier  and  Kezerle,  for  example,  the  gas 

pressure  gradient  term  is  missing  from  the  solid  phase  momentum  equation. 

The  reason  for  this,  as  pointed  out  in  the  next  chapter,  is  that  in  so 

doing  the  equations  become  hyperbolic.  Hughes  [6]  has  shown  that 

3P 

retention  of  the  term  - (1-(J))  ^ in  the  solid  phase  momentum  equation 
results  in  the  mechanical  set  being  non-hyperbolic . The  form  derived 
in  this  study  agrees  with  the  work  of  Kuo  and  Summerfield  [17]  but  not 
with  the  later  work  of  Kuo,  e^  al^. , [22]  where  pressure  was  included 
in  the  term  x^,  nor  with  the  model  of  Gough  [18]  or  Culick  [23]  where 
the  term  is  retained. 

The  pressure  gradient  term  is  also  the  main  source  of  differences 
between  models  in  the  gas-phase  momentum  equation.  All  investigators 
write  this  term  as  the  partial  gradient  of  pressure  except  for  Kuo 
and  Summerfield  [17],  who  write  it  as  the  gradient  of  a partial  pres- 
sure. The  later  model  by  Kuo  e^  aj^. , [22]  subtracts  the  term  - P 
by  redefining  the  drag  coefficient. 

The  final  forms  of  the  equations  given  in  Table  2.1  by  Krier  and 
Van  Tassel  1 [1]  are  distinct  from  the  other  models  in  that  they  were 
derived  from  a continuum  mixture  approach  and  not  the  separated  flow 
approach.  An  additional  term  appears  in  each  momentum  equation  as 
explained  in  Reference  [b]  and  supported  by  Soo  [24].  These  terms 
arise  during  the  splitting  of  the  mixture  momentum  equation  and  are 
called  "inertial  interaction"  terms.  Also  the  typical  source  term 

r (u  - u ) now  appears  in  both  phases  as  T (u 
S S P S 8 

where  u is  a mean  velocity  defined  as 
m 


u ) and  r (u  - u ) 
m-^  g'-  g m^ 


u = 
m 


(Pj  ^ P2 

(Pi  - P2) 


(2.14) 


These  continuum-derived  equations  are  included  here  mainly  for  compar- 
ison. In  general  those  equations  are  more  rigorous. 

The  first  five  sets  of  momentum  equations,  although  derived  by 
apparently  different  means,  have  essentially  arrived  at  the  same  final 
form.  It  is  possible  that  this  means  that  the  averaging  and  control 
volume  techniques  match  the  separated  flow  derivation  as  utilized  by 
Culick  and  as  presented  here  in  Appendix  A. 


9 


The  energy  equations  of  all  six  groups  of  investigators  are  shown 
in  Table  2.2.  The  gas  phase  energy  equations  of  both  Gough  and 
Culick  were  reported  in  convenient  "operator"  (Eulerian)  form  and  were 
copied  directly.  The  equations  of  Kuo  and  his  co-workers  required 
careful  manipulation  to  change  them  into  this  form  from  the  conserva- 
tive form  in  which  they  were  reported. 

Focusing  attention  first  on  the  gas-phase  energy  equations  one 

can  see  that  the  control  volume  technique  utilized  here  yields  equations 

which  are  surprisingly  close  to  those  derived  by  other  approaches,  with 

only  a few  significant  differences.  In  the  model  of  Kuo  and  Summerfield 

there  is  a minor  difference  in  that  the  pressure  work  term  is  written  as 

au  9(<j>u  ) 

- instead  of  - P — . Gough,  Culick,  and  the  later  model  by 

Kuo  use  this  latter,  more  acceptable  form.  .\n  unexplainable  difference 
is  the  additional  term  “ Pg  used  by  Gough  and  Kuo,  Koo,  Davis,  and 
Coates.  Gough  also  has  an  extra  term  Obiter  discrepancies 

between  the  work  done  here  and  that  of  Kuo  e^  al. , arise  from  his 
re-definition  of  the  interphase  drag.  The  result  is  an  additional, 
unusual  term,  u P in  their  equation.  The  model  of  Culick  corre- 
lates  well  with  the  model  used  here,  except  that  he  chooses  to  place 
some  terms  in  the  gas  phase  equation  which  more  properly  belong  in  the 
solid  phase  equation. 

Turning  to  the  solid  phase  energy  balances,  one  sees  immediately 
that  neither  Kuo  and  his  co-workers  nor  Gough  writes  such  an  equation. 
Instead  they  have  resorted  to  calculating  the  particle  surface  temper- 
ature by  solving  the  unsteady  heat  conduction  in  the  solid  phase.  Kuo 
and  Summerfield  have  made  comments  in  defense  of  not  writing  a solid 
phase  energy  equation.  There  is  certainly  some  question  as  to  what 
mixture  energy  is  being  conserved  when  one  does  not  write  such  an 
equation  in  a form  syTiunetrical  to  the  gas  phase  energy  equation  as  done 
here  and  by  Culick. 


2.5  Notes  on  Constitutive  Relations 


In  order  to  solve  the  equations  listed  in  section  2.3  several 
variables  must  be  set  through  the  use  of  known  physical  laws  governing 
the  critical  interaction  processes.  In  particular,  the  following  must 
be  specified. 
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1)  An  ignition  criterion,  in  terms  of  the  bulk  temperatures  of 

the  solid,  T = E /C  . 

P P V 

2) *  The  inter-phase  heat  transfer  coefficient,  h’gp- 

3)  The  inter-phase  drag  coefficient,  f^^. 

4)  The  propellant  burning  rate,  f^. 

5)  The  axial  normal  stress  due  to  particle-particle  interaction, 

T , 


S 

6)  The  variation  of  viscosity  with  temperature. 

The  constitutive  relations  as  used  here,  which  determine  the 
relations  listed  above,  are  given  and  discussed  in  Appendix  B.  It 
can  be  noted  here  that  in  contrast  to  the  earlier  work  of  Dimitstein 
[4]  only  the  stated  drag  laws  were  changed  through  multiplication  by 
a variable  constant  and  not  the  heat  transfer  coefficient.  Also  the 
wall  shear  forces,  heat  losses  to  walls,  and  mass  sources  from  the 
walls  were  not  necessary  and  properly  neglected  in  this  one-dimen- 
sional model. 

Although  little  emphasis  is  placed  on  these  relations  here, 
their  proper  form  and  function  are  critically  important  to  the  success 
of  this  model  for  reactive  two-phase  flow. 


Note  that  'I'  H h (T  -T  ) (p  (S/V)  , where  (S/V)  represents  the 
SP  S P P P 

surface  to  volume  of  the  particle,  which  for  a sphere  = S/r^. 


Here , F 


gP 


(u  - u ) 
g P 


SOLID  PIUSE  COMMENTS 
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CHAPTER  THREE 
NUMERICAL  INTEGRATION 

3. 1 Introduction 

The  two-phase  separated  flow  model  generates  a set  of  six  non- 
linear, inhomogeneous,  and  coupled  hyperbolic  partial  differential 
equations,  as  demonstrated  in  the  previous  chapter.  The  numerical 
method  chosen  for  solving  these  equations  is  the  explicit  two-step 
MacCormack  scheme.  This  method  was  chosen  for  its  accuracy  and  its 
ability  to  capture  shocks  without  the  use  of  "artificial  viscosity." 

With  an  explicit  numerical  method  it  is  possible,  under  certain 
circumstances,  to  "march  out"  a solution  to  a system  of  hyperbolic 
equations  on  a fixed,  rectangular  grid  instead  of  a grid  composed 
of  the  characteristics  of  the  equations.  For  the  "two-dimensional" 
system  of  concern  here  (one  dimension  in  space  and  one  in  time)  a 
grid  of  the  following  type  is  constructed. 


Figure  3.1 

As  discussed  in  Ames  [25],  the  solution  to  the  hyperbolic  system  is 
known  within  the  triangle  formed  by  the  real  characteristics,  lines 
AC  and  BC,  in  Fig.  3.1  if  the  solution  is  known  up  to  the  time  indi- 
cated by  the  line  AB.  The  point  predicted  by  the  numerical  computa- 
tion of  one  step  forward  in  time,  point  D,  is  a product  of  the  grid 
spacings  which  determine  the  slope  of  the  "finite  difference  character- 
istics." If  the  slope  of  these  finite  difference  characteristics  is 
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greater  than  the  slope  of  the  real  characteristics,  instability  sets 
in;  if  equal,  the  solution  is  identical  to  the  characteristic  solu- 
tions; if  less  than,  the  system  is  stable  but  accuracy  decreases  with 
decreasing  slope.  The  idea  of  an  explicit  finite  difference  method, 
then,  is  to  accurately  approximate  the  real  characteristics  of  the 
system  on  a rectangular  grid.  These  points  must  be  considered  later 
when  the  question  of  stability  is  studied. 

3.2  The  MacCormack  Scheme 

The  MacCormack  scheme  is  a two-step,  predictor-corrector  method 
of  integration.  The  basic  idea  of  the  method,  as  applied  to  partial 
differential  equations  can  readily  be  understood  through  an  analogy 
with  an  ordinary  differential  equation.  Consider  the  simple  O.D.E. 

It  = • (3 

Knowing  the  value  of  the  function  u at  a time  n,  the  value  of  a new 
time,  n + 1,  is  sought. 


I L. 

n n+1  t 

Figure  3.2 

An  approximate  value  at  the  new  time  is  found  in  the  predictor 
step,  knowing  only  the  value  of  the  function  and  its  slope  at  time  n. 

r-  rj 

n+I  n du  ,, 

“ = ” ^ ^ndTj  • f- 

The  bar  above  a term  is  used  to  indicate  the  intermediate  nature  of 
a value  computed  in  the  predictor  step.  As  shown  in  Fig.  3.3,  this 
value  can  be  in  considerable  error. 
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n n+1  t 


Figure  3.3 

The  "predicted"  value  is  revised  by  first  finding  a new  slope  of  the 
function  at  the  predicted  point.  The  "corrected"  value  of  the  func- 
tion at  time  n + 1 is  then  found  by  averaging  the  known  slope  of  the 
function  at  time  n and  the  predicted  slope  of  time  n+1  and  making 
a new  extrapolation  from  the  known  value,  u^. 


n n+1  t 

Figure  3.4 

The  equation  for  finding  this  more  accurate  point  is: 


If  a substitution  for  the  slope  of  the  function  at  time  n + 1 is 
made  in  Eq.  (3.3)  from  Eq.  (3.2),  another  common  form  of  the  corrector  1 

step  is  obtained.  ’ 

(3.4  I 


Figure  3.5 

The  application  of  the  MacCormack  scheme  to  a partial  differen- 
tial equation  is  a natural  extension  of  the  method  shown  above.  In 
this  case,  the  simple  vector  form  of  the  equation  is; 


3u 

at 


+ 


0 


in  - a If. 
at  ■ ^ ax 


(3.5) 


(3.6) 


where  u and  F are  given  column  vectors  representing  the  parameters  of 
interest  (p,  u,  and  e)  and  A(=  - B)  is  a constant  matrix  of  coefficients 
The  formulation  of  the  predictor  and  corrector  steps  is  the  same  as  be- 
fore, but  now  the  time  derivatives  of  the  function  u can  be  replaced  by 
special  derivatives  in  the  coordinate  x as  Eq.  (3.6)  clearly  shows.  In 
the  simplified  case  u = F,  the  following  grid  system  can  be  constructed. 
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Initial  Conditions 
Figure  3.6 


The  equation  for  the  predictor  step  is  now  a revision  of  Eq.  (3.2). 
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(3.7) 


A forward  difference  in  space  was  chosen  to  evaluate  the  spatial 
derivative,  although  it  was  not  mandatory.  It  follows  that  the  equa- 
tion of  the  corrector  step  is  a simple  revision  of  either  Eq.  (3.3) 
or  Eq.  (3.4).  Choosing  Eq.  (3.4)  yields; 
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When  both  steps  are  completed,  the  value  of  the  function  u is 
known  for  a new  time  at  one  location.  (See  Fig.  3.6.)  The  entire 
process  of  both  predictor  and  corrector  is  then  repeated  for  the  next 
point  in  the  x-direction  and  the  same  time.  When  all  the  points  in 
the  one-dimensional  grid  are  known  at  the  new  time  another  step  forward 
in  time  can  be  taken  and  the  entire  process  is  repeated  for  each  point 
starting  at  a known  boundary. 

Although  each  step  (predictor  and  corrector)  of  the  process  is  a 
simple  forward  or  backward  finite  difference  of  first  order  accuracy, 
the  overall  process  has  second  order  accuracy  as  explained  by  Warming 
[26].  This  procedure  also  acts  as  an  implicit  artificial  viscosity  to 
help  capture  embedded  shocks.  The  above  completes  description  of  the 
scheme  and  reasons  for  its  choice. 

3.3  Choice  of  Mesh  System  and  Boundary  Conditions 

There  are  two  basic  mesh  systems  to  which  the  MacCormack  scheme 
can  easily  be  applied.  (See  Fig.  3.7  below). 


Figure  3.7 

The  major  difference  between  the  two  systems  is  that,  in  the  first, 
there  is  a grid  point  lying  on  the  wall,  while  in  the  second  there  is 
no  such  point.  The  choice  of  a mesh  system  has  a great  bearing  on  how 
boundary  conditions  are  employed.  For  the  work  done  here  a mesh  sys- 
tem of  the  first  type  was  used  and  labeled  as  shown  in  Fig.  3.8.  This 
convention  must  be  kept  in  mind  as  the  boundary  conditions  are  discussed. 
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Figure  3.8 

Since  both  ends  of  the  one-dimensional  flow  were  considered 
closed  a solid  wall  condition  was  imposed  at  both  places.  The  only 
boundary  condition  that  can  be  stated  with  certainty  is  that  the 
velocity  of  both  particles  and  gas  must  be  zero  at  the  end  walls 
[27].  The  most  convenient  way  of  applying  this  condition  is  to  de- 
fine a row  of  fictitious  points  inside  the  wall  and  use  an  anti-sym- 
metric  reflection  of  velocity,  u,  while  also  setting  u = 0 at  the 
walls.  The  procedure  is  very  compatible  with  the  first  mesh  system. 
Boundary  conditions  on  the  other  parameters  are  not  as  easy  to  deduce. 

In  order  to  derive  a boundary  condition  on  gas  temperature  an 
assumption  must  be  made  as  to  the  nature  of  the  wall.  One  choice  is 
to  assume  an  adiabatic  wall,  although  that  was  not  done  here.  Instead 
it  was  assumed  that  the  heat  conducted  by  the  wall  is  equal  to  the  heat 
convected  to  the  wall  by  the  gas  as  was  done  by  Krier,  Dimitstein,  and 
Gokhale  [5].  Thus 

9T 

h (T  - T ) = K ^ , where  the  heat  transfer  coefficient,  (3.11) 
^ ^ g ^ h^,  is  an  assumed  given  function. 

At  the  right-hand  wall,  for  example,  this  becomes 


h (T 
c g 


- T ) 
w 


■ " '^N2^ 


(3.12) 


This  can  be  solved  for  T (=T.,,). 

V N1 


Since  the  Noble-Abel  non-ideal  equation  of  state  is  assumed  to 
hold  at  all  points,  the  remaining  boundary  condition  on  the  gas  phase 
need  only  specify  either  pressure  or  density  and  the  other  is  de- 
termined. Because  this  state  equation  employs  a variable  co-volume 
which  is  a function  of  density  it  was  necessary  to  impose  a boundary 
condition  on  density  so  that  with  density  known,  the  co-volume  could 
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be  derived,  and  then  the  pressure  calculated.  The  condition  imposed 
at  the  wall  was  = 0 as  suggested  by  Dwyer  and  Allen  [28].  This 

meant,  at  the  left-hand  wall  for  example: 

p(l)  = p(2)  = p(3)  . (3.13) 

The  boundary  values  for  the  solid  phase  were  calculated  using 
anti-synunetric  reflection  of  velocity  (as  stated  earlier),  symmetric 
reflection  of  other  variables,  and  employing  regular  interior  point 
differencing. 

To  summarize,  the  boundaries  were  handled  in  the  following  way: 


Clas  Phase: 


= 0 


= 0 


= f (heat  balance) 


(3.14) 


Solid  Phase: 


P2(l)  = P2(3) 


U2(l)  = - U2(3) 


e2(l)  = e2(3) 


P2(N)  = P2(N2) 


U2(N)  = - U2(N2) 


e2(N)  = e2(N2) 


(3.15) 


where  p^  = <J)Pg  and  p2  = (l-<J))Pp  as  defined  earlier,  ^lnd  where  [ ® 

gas  and  [ ]2  “ solid  phase. 

Values  for  the  gas  phase  variables  at  the  fictitious  points  were 
also  assigned  through  reflection.  It  must  be  noted  however  that  these 
values  were  employed  only  in  the  calculation  of  the  MacCormack  damping 
factor,  as  discussed  later.  They  were  not  used  to  find  gas  phase 
boundary  values.  There  has  been  some  concern  expressed  [29]  that  the 
reflection  principle  implies  "extra"  boundary  conditions. 


9f 

57 


= 0 


p.u,  c 


(5.16) 
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It  is  felt  that  this  is  not  true  for  the  equations  studied  here,  due 
to  the  large  source  and  sink  terms.  It  must  also  be  noted  that  with 
the  mesh  system  used  reflection  through  a wall  does  not  insure  zero 
gradient  at  the  wall. 

3.4  Initial  Conditions 

Initial  conditions  specified  temperature,  pressure,  and  velocity 
throughout  the  grid.  Temperatures  were  set  by  a linear  decay  over  the 
first  five  grid  points.  Pressure  was  subject  to  an  exponential  decay 
throughout  the  entire  domain.  Velocities  of  gas  and  solid  at  time  = 0 
were  set  equal  to  zero  everywhere.  Initial  porosity  was  set  as  de- 
sired, usually  distributed  evenly.  These  conditions  are  shown  graphi- 
cally in  Fig.  3.9. 


Figure  3.9 
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The  pressure  profile  initialized  above  differs  from  that  used  by 
Dimitstein  f4J  where  the  exponential  decay  was  applied  to  only  one 
quarter  of  the  grid  points.  It  was  found  in  the  course  of  this  work 
that  a "smooth"  initial  pressure  profile  (no  discontinuities)  was 
needed  to  provide  stable  results. 


3.5  Discussions  Regarding  Numerical  Stability 

The  ability  of  an  explicit  finite  difference  method  is  dependent 
upon  the  grid  size  and  time  step  as  mentioned  in  Section  3.1.  It  has 
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been  shown  that  in  order  to  obtain  a stable  solution  to  the  hyperbolic 
system  the  numerically  predicted  point  must  lie  within  the  area  formed 
by  the  real  characteristics  (See  Fig.  3.1).  For  the  MacCormack  scheme 
this  is  assured  by  satisfying  the  Courant,  Friedrichs,  Levy  criteria 
[30] , i.e. , 


At(c  -t-  |u|)  < 
Ax 


(3.17) 


where  c represents  the  average  speed  of  sound  of  the  system. 

With  a fixed  grid  size,  a slightly  more  restrictive  form  of  this 
criteria  was  used. 


.7*Ax 
(c  + lul) 


(3.18) 


There  still  exists  a certain  amount  of  ambiguity  when  applying 
this  restriction  to  a two-phase  mixture.  The  question  arises  whether 
the  velocities  c and  |u|  should  refer  to  the  mixture  or  just  to  the 
gas.  This  dilemma  is  crucial  since  it  is  well  known  that  the  speed  of 
sound  in  a solid-gas  two-phase  mixture  can  be  less  than  the  speed  of 
sound  in  either  phase  alone.  For  this  research,  stable  results  were 
obtained  when  c was  chosen  as  the  mixture  sound  speed  and  |u|  was  taken 
as  the  highest  gas  velocity  encountered.  An  expression  for  the  sound 
speed  of  a two-phase  mixture  can  be  found  in  Wallis  [8] , as 


1 f. 

-2=  ap  + (1-a)  ^ 

L Jl2  2 ^1  1_ 


(3.19) 


It  can  be  seen  that  the  sound  speed  of  the  mixture  is  a function  of 
the  sound  speed  in  each  phase  separately.  Therefore  the  first  priority 
was  to  find  a value  for  the  sound  speed  in  each  phase.  For  the  gas 
phase,  standard  isentropic  sound  speed  was  used,  namely 


= /yOT 


(3.20) 


where  T was  chosen  as  the  initial  gas  temperature.  For  the  solid 
phase , 
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where  is  the  bulk  modulus  of  a medium  hard  plastic.  These  calcula- 
tions are  carried  out  in  Appendix  C.  With  these  values  known,  the 
sound  speed  of  the  mixture  was  calculated  from  Eq.  (3.16).  The  value 
computed  was  c = 2,708  ft/sec, which  was  less  than  the  speed  of  sound 
in  either  phase  alone. 

With  c and  |u|  determined,  a desirable  time  step  was  calculated 
from  the  C.F.L.  criteria  Eq.  (3.15).  It  was  found  that,  using  this 
restriction,  stable  solutions  were  always  obtained  for  those  systems 
of  equations  believed  to  be  hyperbolic.  A recent  paper  by  Krier, 
Gokhale,  and  Hughes  [6]  discusses  the  concern  and  manner  in  which  the 
governing  equations  are  actually  hyperbolic. 
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CHAPTER  FOUR 
THE  PREDICTED  RESULTS 

4. 1 Introduction 

As  discussed  in  the  previous  chapters,  the  major  purpose  of  this 
work  is  to  predict  the  pressure  wave  fronts  relative  to  the  ignition 
front  in  a two-phase  reactive  mixture  in  order  to  give  some  indication 
whether  a high-speed  transient  deflagration  wave  might  transit  to  a de- 
tonation in  a highly  loaded  propellant  bed.  To  properly  study  all  pos- 
sible variations  in  the  allowed  parameters,  the  scope  of  this  work  would 
have  had  to  be  broadened  beyond  its  original  intent.  This  not  being 
possible  in  the  time  period  available,  only  a limited  number  of  results 
are  included  here.  It  is  expected  that  this  work  will  be  continued  and 
broadened  in  the  near  future. 

After  several  hundred  independent  calculations  a baseline  case  with 
standard  input  conditions  was  deduced.  The  governing  equations  solved 
are  those  shown  in  Chapter  Two,  namely  Eqs.  C2.1)  to  (2.7),  with  con- 
stitutive relations  as  outlined  in  Appendix  B.  Specifically,  the  solid 
phase  stress  tensor  in  Eq.  (2.4)  does  not  include  the  gas  pressure.  As 
mentioned  previously,  the  reason  for  this  is  that  it  has  been  shown  [6] 
that  the  set  of  equations  will  be  totally  hyperbolic  only  if  there  is  no 
gas  pressure  gradient  in  the  solid  phase  momentum  equation. 

Table  4.1  summarizes  the  most  important  input  used  to  carry  out 
the  typical  results  shown  below.  The  main  results  presented  from  this 
analysis  are  the  pressure  distributions  as  time  progresses,  the  flame 
front  locus  and  pressure  front(s),  as  well  as  developments  in  the  tem- 
perature and  velocity  fields. 

4. 2 Results  From  the  Baseline  Case 

Figures  4.1  through  4.5  present  the  distribution  histories  of  the 
fluid  dynamic  variables  in  a bed  of  solid  propellant,  3 inches  long, 
with  a solids  loading  at  60  percent  ((j)^  = 0.40).  The  results  in  these 
figures  were  termed  the  baseline  case,  as  defined  above. 

Figure  4.1  outlines  the  calculated  pressure  distributions  at  four 
different  times.  As  in  previous  work  related  to  this  problem  (see 
Krier  [3,5],  the  appearance  of  a "continental-divide"  in  the 
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Table  4.1 


TYPICAL  INPUT  DATA 


Paramet er 

Value 

Initial  bed  temperature,  T^,  T^,  T^ 

SSO^R 

Ignition  temperature, 

545“r 

Initial  bed  porosity,  (f) 

0.40 

Propellant  burning  rate  constant,  bj 

0.0  in/ sec 

Propellant  burning  rate  proportionality  constant,  b^ 

0.00553  in/sec-Psia*' 

Propellant  burning  rate  index,  n 

(0.90) 

Propellant  density, 

0.0571  Lbm/in^ 

Initial  grain  diameter,  d 

400ym  (200  ym) 

Chemical  energy  released, 

2360.9  BTU/Lbm 

Molecular  weight  of  gas,  MW 

22.6  Lbm/Lbmole 

Covolume  of  propellant  gas, 

29.85  inVhbm  . | 

Specific  heat  ratio  of  gas,  y 

1.252 

Gas  viscosity,  yg 

0.249  X 10"^lbm/in-sec 

Universal  gas  constant,  R 

1.9869  BTU/Lbmole-‘’R  1 

Thermal  conductivity  of  the  walls, 

0.3  X 10"^  BTU/in-sec-®R  U 

•( 

Total  bed  length,  Lg 

3.0  in.  j| 

Percent  heat  transfer  to  particle  and 

10%  Ij 

to  wall  after  ignition 
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interior  of  the  bed  is  pronounced  and  pressure  gradients  are  ex- 
tremely steep.  In  addition,  the  average  pressures  behind  the  front 
are  very  high,  reaching  in  excess  of  O.S  x 10°  psia.  It  appears 
that  for  nearly  the  same  constituitive  laws,  the  separated  flow  model 
(as  derived  here)  predicts  pressure  wave  mechanics  similar  to  that 
calculated  with  the  continuum  mixture  model  of  Krier  et  al . [6] . 

The  locus  of  the  flame  (ignition)  front  for  the  baseline  case 
is  presented  in  Figure  4.2.  As  expected,  the  deflagration  speed 
(slope  of  the  x-t  locus)  accelerates  to  several  mm/y-sec  in  this 
3-inch  long  bed.  Also  shown  in  this  figure  is  a line  labeled  the 
"pressure  front."  This  locus  was  derived  by  noting  the  midpoint  of 
the  pressure  wave  shown  in  Figure  4.1  at  various  times.  It  appears 
that  the  pressure  front  and  the  flame  front  are  coincident  as  the 
flame  accelerates  to  complete  ignition. 

Figure  4.3  shows  both  the  gas  velocity  and  the  propellant  par- 
ticle velocity  distributions  at  the  chosen  times.  As  might  be  ex- 
pected, the  velocity  peaks  shown  are  a consequence  of  the  steep  pres- 
sure front  arriving  at  the  various  locations.  Also  worth  noting  are 
the  extreme  gas  velocities,  often  exceeding  5000  ft/sec,  and  the  fact 
that  peaks  in  particle  velocity  lag  peaks  in  gas  velocity  at  any  given 
time. 

Figure  4.4  presents  gas  temperature  and  particle  temperature 
distribution  histories.  With  the  chemical  energy  fixed  for  this 
propellant,  the  gas  temperatures  are  generated  at  values  of  roughly 
6500®R.  These  temperatures  decay  as  the  gases  lose  heat  to  the  solid 
by  convection  into  the  bed  interior.  Later,  as  the  pressure  front 
steepens,  the  gases  are  compressed  in  the  interior  to  values  far  ex- 
ceeding their  incoming  values.  The  subsequent  particle  temperatures 
are  increased  not  only  due  to  the  convective  heat  transfer,  but  also 
due  to  the  rapid  changes  in  kinetic  energy  of  these  mobile  particles. 

The  last  predictions  displayed  for  the  baseline  case  are  the 
porosity  distribution  histories  shown  in  Figure  4.5.  The  porosity 
profiles  show  that  near  the  front  of  the  bed  (which  is  partially  ig- 
nited at  t = 0)  the  porosity  increases  fairly  rapidly.  This  is  due  to 
both  a reduction  in  volume  of  the  particles  due  to  burning  and  the 
forward  motion  of  these  particles  being  dragged  along  by  the  gases. 
Subsequent  to  this,  the  particles  are  compacted  somewhere  in  the  bed 


Pressure  ( Kpsi ) 


Figure  4.1  Pressure  distribution  history  predicted  for  a three  inch 
long  bed  with  an  initial  solids  loading  of  60  percent. 
[Predictions  for  Figs.  4.1  - 4.5  utilized  the  second-form 
of  the  particle  energy  equation.] 
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Figure  4.3  Gas  velocity'  and  particle  velocity  distribution 
history. 
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Figure  4.4  Particle  temperature  and  gas  temperature  distribution  history. 
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Figure  4.5  Porosity  (gas  volume  fraction)  distribution  history 
complementary  to  predictions  in  Figs.  4.1  - 4.4. 
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interior  to  solids  loadings  greater  than  the  original  60  percent. 
(Recall  that  the  solid  fraction  equals  !-<).)  For  the  results  shown 
here,  a significant  fraction  of  the  bed  has  been  consumed  at  a time 
equal  to  50  y-sec. 


4.3  Parametric  Variations  on  the  Baseline  Case 

As  discussed  in  Section  4.1,  it  was  not  possible  to  vary  all  input 
parameters  or  study  all  possible  formulations  of  the  governing  equations. 
The  following  figures,  however,  will  indicate  the  sensitivity  of  the 
predictions  shown  in  Figs.  4.1  - 4.5  to  some  of  the  most  important  pa- 
rameters. Although  one  might  generally  assume  that  the  gas  phase  fluid 
and  thermodynamic  variables  are  well  known,  it  was  surprising  to  see  the 
significant  changes  in  the  pressure  wave  history  when  different  relations 
are  used  to  determine  gas  viscosity.  Figure  4.6  clearly  indicates  these 
differences  for  the  two  relations  assumed  to  represent  gas  viscosity  as 
a function  of  temperature.  It  must  be  remembered  that  viscosity  is  an 
important  parameter  in  both  the  heat  transfer  and  viscous  drag  correla- 
tions used  here.  The  baseline  case  (as  shown  as  a dotted  curve)  uses 

y = 2.49  X 10'®  (T/T^)-®®  Ibm/in-sec  . (4.1) 

When  an  alternate  relation,  namely. 


7.57  X 10  **  (T/Tq)^-® 
(T  + 198) 


C4.2) 


was  used,  these  interphase  correlations  changed  so  drastically  that  re- 
latively long  delays  resulted  before  the  steep  pressure  fronts  de- 
veloped. Two  profiles  (shown  as  solid  lines)  are  presented  in  Fig.  4.6 
for  this  second  relation.  This  delay  is  also  apparent  in  the  flame  front 
history  as  shown  for  both  cases  in  Fig.  4.7 

Increasing  the  solids  loading  results  in  more  rapid  wave  dynamics 
and  greater  peak  pressures,  as  expected.  Figure  4.8  presents  the 
flame  front  locus  for  three  solids  loadings  varying  from  70  percent  to 
50  percent.  Had  calculations  carried  past  a bed  length  of  three  inches, 
possible  transitions  to  detonation  might  have  been  expected  in  much 
shorter  run-up  lengths  for  the  high  solids  loading.  Shown  as  an  in- 
sert on  Fig.  4.8  is  the  peak  pressure  in  the  bed  interior  at  times  just 
before  full  bed  ignition.  Note  that  increasing  the  solids  loading 


34 


Location  ( inches) 
X — - 


Figure  4.6  Comparison  of  the  pressure  distribution  for  an 
alternate  gas  viscosity  relation  with  that  pre- 
dicted at  an  earlier  time  for  the  standard  viscosity 
law  (dashed  line) . 
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from  50  percent  to  70  percent  triples  this  peak  pressure. 

Decreasing  the  particle  size  from  the  standard  200  y diameter 
(r^  = .004  inches)  yields  a trend  similar  to  increasing  the  solids 
loading.  In  the  case  shown  in  Figure  4.9,  decreasing  the  radius  to 
100  y results  in  an  increase  in  the  S/V  ratio  of  the  particles  at  100 
percent.  (Since  (S/V)^  = 3Cl-6^)/rp)  For  the  input  data  base  used 
here,  the  predicted  pressures  for  the  small  particle  size  exceeded 
10^  psi  at  which  time  calculations  were  terminated.  Therefore  the 
maximum  pressure  in  the  bed  as  a function  of  particle  size,  as  pre- 
sented in  the  insert,  is  shown  for  an  arbitrary  time  of  26.6  ysec, 
and  not  when  the  whole  bed  was  ignited. 

4.4  Results  Due  to  Changes  in  the  Governing  Equations 

There  are  three  current  controversies  that  have  appeared  in  the 
literature  by  investigators  having  developed  approximate  separated 
flow  models.  The  first  deals  with  the  appropriate  form  of  expressing 
the  pressure  gradient  in  the  gas-phase  momentum  equation.  Some  have 
used  the  concept  that  the  term  should  be  written  as  the  gradient  of 
partial  pressure,  i.e.,  - 8/ 9x((l)Pg) . Recall  that  in  the  derivation 
used  here  the  gas  phase  was  subject  to  its  partial  gradient  of  pres- 
sure, -(f)*  3/3x(Pg).  Table  2.1  indicates  what  form  the  other  models  have 
taken. 

To  determine  the  effects  of  using  the  alternate  form  calculations 
were  carried  out  in  which  the  term  -Pg0(})/3x)  was  added  to  the  right 
hand  side  of  the  gas  phase  momentum  equation.  Figure  4.10  clearly 
shows  that  for  the  data  base  used  here  such  an  addition  has  a very 
dominant  effect  on  the  pressure  profile.  One  distribution  (at  t = 42.4 
y-sec)  is  shown  and  compared  with  profiles  that  were  predicted  without 
the  additional  term.  The  flame  front  locus,  however,  is  not  markedly 
different  from  the  standard  result  shown  in  Figure  4.2.  Figure  4.11 
presents  the  flame  front  and  compares  it  with  the  press'ire  front  cal- 
culations at  four  discrete  locations.  The  onset  of  a very  rapid  pres- 
sure rise  was  defined  as  the  arrival  of  this  pressure  front. 

The  second  possible  variation  in  the  governing  equations  concerns 
the  possible  inclusion  of  the  partial  gradient  of  gas  pressures  in  the 
solid  phase  momentum  equation,  namely  - (l-(ti)3  Pg/3  x.  Hughes  [6]  has 
shown  that  this  term  in  that  equation  results  in  a non-hyperbolic 
system  and  one  that  might  experience  difficulties  in  numerical 


The  flame  front  locus  predicted  for  the  case  when  the  term,  - d/dx(  <PP) , is 
used  in  the  gas-phase  momentum  balance. 
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integration.  As  shown  in  Table  3.1,  several  modelers  retain  such  a 
term.  The  calculations  to  date  seem  to  confirm  the  concepts  of 
Hughes,  although  the  evidence  is  not  conclusive.  Stable  and  well- 
behaved  results  were  never  really  achieved  when  this  term  was  in- 
cluded. Pressure  and  velocity  oscillations  usually  occurred  similar 
to  the  problem  that  plagued  Dimitstein  and  Krier.  Altering  the  mesh 
size  smoothed  results  somewhat.  Use  of  a damping  factor  recently  de- 
vised by  Hung  and  MacCormack  [31]  had  little  effect.  Settling  this 
controversy  will  require  many  more  additional  calculations. 

A third  variation  on  the  separated  flow  model  is  concerned  with 
the  formulation  of  the  solid  phase  energy  equation.  Appendix  A pre- 
sents discussion  on  why  it  is  possible  to  have  three  separate  forms. 

In  addition,  some  investigators  do  not  bother  to  solve  any  such 
equation  (see  Table  3.2).  Figure  4.12  shows  the  pressure  distribution 
histories  for  the  three  energy  equations  derived  in  Appendix  A.  The 
first  two  forms,  labeled  (a)  and  (b) , differ  only  in  the  manner  in 
which  the  pressure-work  term  is  expressed.  Recall  that  the  baseline 
case  utilizes  the  second  form,  (b) . The  third  form,  (c)  does  not 
assume  that  the  particles  are  fluid- like.  That  is,  there  is  no  pos- 
sible interchange  between  kinetic  and  internal  energy.  These  pressure 
profiles  were  calculated  with  the  second  viscosity  law.  It  is  ap- 
parent from  these  results  that  much  more  work  will  be  required  to 
justify  which  of  these  equations  should  be  used.  At  this  point,  the 
first  form  of  the  solid-phase  energy  equation  was  discarded  because 
predicted  ;:Atrosity  profiles  showed  severe  compaction  and  gas  tempera- 
tures were  much  too  high. 

Figures  4.13  through  4.16  present  results  using  the  third  form  of 
the  energy  equation,  namely  equation  16  in  Appendix  A.  Figure  4.13 
presents  the  pressure-time  history  at  six  locations,  clearly  showing 
the  development  of  a pressure  front.  The  gas  temperature  and  particle 
temperature  distribution  histories  are  shown  in  Figure  4.14  and  4.15. 
The  calculation  of  most  interest  is  shown  in  Figure  4.16.  Notice  the 
abrupt  change  in  flame  speed  at  approximately  60  p-sec  which  coincides 
with  the  intersection  of  the  pressure  front  obtained  from  Figure  4.13. 

The  results  shown  in  Figure  4.16  might  represent  a possible  de- 
flagration to  detonation  transition.  Beyond  60  y-sec  the  pressure 
front  leads  the  flame  front  as  one  could  expect  to  see  in  a detonation 
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Figure  4.12  Pressure  distribution  at  a given  time  for  three  different 
forms  of  the  particle-energy  balance.  [The  alternate  gas 
viscosity  law  was  utilized  for  all  three  calculations.] 


Figure  4.13  The  pressure-tirae  variations  at  six  different  locations. 
jFigs.  4.13  - 4.16  utilize  the  third  form  of  the 
particle-energy  equation.] 


2.5 


3. 


Location  (inches  ) 


Figure  4.14  Gas  temperature  distribution  history. 
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Figure  4,16  Flame  front  locus  and  pressure  wave  front,  indicating  conditions 
for  possible  DDT. 
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shock  supported  by  a combustion  region.  Also,  the  flame  front 
accelerates  from  a speed  of  1.9  mm/ysec  at  x = 2 inches  to  about 

2.5  nun/y-sec  at  x = 2.75  inches.  It  is  concluded  that  a search  for 
possible  DDT  should  include  a thorough  study  of  this  alternate  form 
of  the  solid-phase  energy  equation. 

Other  details  are  shovm  in  Fig.  4.17  again  depicting  the 
pressure  distribution  history  for  a modified  baseline  case.  Here  the 
drag  is  twice  that  used  in  all  previous  calculations.  This  figure  is 
included  to  indicate  that  a steep  pressure  reversal  is  predicted  after 
the  pressure  front  strikes  the  rear  wall.  These  large  reverse  gra- 
dients are  expected.  It  is  possible  that  additional  gas  compression 
would  result  in  an  alternate  reason  for  possible  DDT.  More  work  is 
required  to  clarify  this.  It  can  be  concluded,  however,  that  the 
MacCormack  scheme  is  quite  stable  for  this  model  and  results  indicate 
pressure  fronts  moving  until  all  particles  are  consumed.  It  is  be- 
cause of  this  result  that  one  can  be  confident  that  the  integration 
scheme  used  here  produces  accurate  and  stable  results  for  the  sepa- 
rated flow  model. 

4.5  Other  Sensitivity  Studies 

Figure  4.18  presents  the  locus  of  flame  front  for  various  ignition 
energies,  or  bulk  ignition  temperatures.  As  expected,  the  lower  bulk 
ignition  temperature  results  in  a shorter  time  required  for  complete  bed 
ignition.  However,  the  ignition  temperature  does  not  seem  to  affect  the 
flame  speed,  i.e.,  the  slope,  on  this  x-t  diagram.  It  is  also  interesting 
to  note  the  higher  maximum  gas  pressure  in  the  bed  interior  for  higher 
ignition  temperatures.  (See  insert  in  Fig.  4.18.)  A simple  explanation 
could  be  as  follows.  The  time  required  for  igniting  an  equal  length  of 
bed  will  be  greater  for  a higher  ignition  temperature,  since  more  energy 
needs  to  be  transferred  to  the  particles  from  the  hot  gases.  In  this 
extra  time  required,  the  portion  of  the  bed  already  ignited  keeps  on 
generating  more  hot  gases.  This  gas  production  becomes  especially  sig- 
nificant for  propellants  with  a high  burning  rate  index  (n  = 0.9).  We 
can  conclude  that,  at  the  conditions  being  studied  here,  the  ignition 
energy  (or  temperature)  variation  is  not  necessarily  of  primary  im- 
portance for  setting  up  transition  to  detonation.  Of  course,  for 
quantitative  analysis,  such  as  a detonation  run-up  length  the  ignition 
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Figure  4.17  Pressure  distribution  history,  indicating  strong  pressure 
wave  reversal  after  the  pressure  front  reflects  from  rear 
wall.  [Viscous  inter-phase  drag  assumed  here  is  twice 
that  used  for  calculations  shown  in  Fig.  4.1.] 
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temperature  (energy)  should  be  specified  as  accurately  as  possible. 

Figure  4.19  depicts  the  effect  of  drag  correlation  on  the  flame 
locus.  In  general,  the  drag  coefficient,  modifying  the  function 
f^p  has  minimum  effect  on  either  the  time  required  for  complete  bed 
ignition  or  the  flame  speeds.  But  it  is  interesting  to  note  that  the 
full  Ergun's  drag  (C^  = 1)  did  not  compute  for  this  separated  flow  model. 
It  was  observed  that,  even  with  = 5.0,  only  50%  of  the  bed  was  ig- 

nited before  the  flow  gradients,  in  particular  the  pressure  gradient, 
became  unrealistically  large,  and  the  program  became  unstable.  A 
larger  value  of  means  a lower  drag.  For  lower  drag  conditions, 
the  particles  are  not  accelerated  enough  to  cause  localized  compaction, 
and  hence  the  maximum  pressures  predicted  were  significantly  lower.  Or, 
corresponding  to  lower  drag,  the  gases  will  travel  faster  and  convective 
heat  transfer  is  more  rapid  with  the  pressure  more  uniform. 

Figure  4.20  shows  the  pressure  distribution  in  a three-inch  bed  at 

one  typical  instant.  Here  H is  defined  as  the  ratio  of  actual  gas-solid 

heat  transfer  coefficient  ^o  that  calculated  by  the  Denton  formula.  The 

lower  the  value  of  H,  the  less  heat  will  be  transferred  from  hot  gases 

to  the  particles.  That  is,  for  smaller  values  of  H,  less  energy  will  be 

removed  from  the  gas  phase  and  transferred  to  the  particle  phase.  As  a 

result  of  retention  of  energy  content  in  the  gas  phase,  the  pressure 

predicted  at  the  equal  time  will  be  higher.  It  seems  apparent  that  both 

the  interphase  viscous  drag  correlation,  f^^,  and  the  interphase  heat 

transfer  correlation,  h , are  sensitive  functions  of  the  dynamic  flow 

Pg 

being  modelled  here.  There  is  no  data  available  that  indicates  whether 
any  of  the  formulas  used  here,  generally  developed  from  cold-flow, 
steady,  low-pressure  experiments  are  applicable  to  the  transient  high 
pressure  conditions  being  studied  here. 

Figure  4.21  depicts  the  flame  front  (ignition)  locus  again  for  three 
different  heat  transfer  coefficients,  H.  Note  that  the  greater  the  gas- 
particle  heat  transfer,  the  greater  the  time  required  for  the  complete 
bed  ignition.  It  seems  an  extension  of  the  previous  explanation.  Higher 
H means  more  energy  will  be  transferred  to  the  particle  phase  from  the 
gas  phase.  Thus,  the  gas  phase  energy  available  for  convecting  towards 
the  unignited  particles  will  be  less  and,  consequently,  the  time  re- 
quired for  the  complete  bed  ignition  will  be  more.  However,  there  does 
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not  seem  to  be  any  change  in  relative  flame  acceleration  for  various 
heat  transfer  coefficients.  Additionally,  the  predicted  maximum  pre- 
sure  (see  insert)  for  various  heat  transfer  coefficients  fails  to  provide 
insight  into  the  dynamic  behavior  during  the  convection  process. 

Figure  4.22  shows  the  flame  front  locus  for  alternative  forms  of 
the  propellant  burning  rate  law.  An  attempt  was  made  to  make  the  rate 
sensitive  to  the  particle  temperature,  as  shown  by  the  equation  given 
in  the  figure.  Even  though  the  temperature  sensitivity  index,  assumed 
here  to  be  1/4,  is  somewhat  arbitrary,  it  was  felt  that  some  form  of 
temperature  sensitivity  of  rate  law  should  be  accounted  for,  since  pro- 
pellants are  known  to  burn  more  rapidly  as  their  mean  (storage)  tem- 
perature is  increased.  This  is  a first  attempt  in  that  direction.  As 
expected,  the  more  sensitive  the  rate  law  is  to  the  particle  temperature, 
the  faster  the  whole  bed  is  ignited.  This  seems  logical  since  a more 
temperature-sensitive  burning  rate  law  generates  more  hot  gases  which 
will  be  convected  forward  to  produce  higher  gradients.  This  also  ex- 
plains the  higher  maximum  bed  interior  pressure  for  the  more  tempera- 
ture-sensitive burning  rate  law,  as  shown  in  the  insert  of  the  figure. 

Finally,  Fig.  4.23  shows  the  flame  front  loci  for  various  coef- 
ficients multiplying  the  particle-particle  interaction  law.  In  general, 
it  shows  that  the  lower  the  particle-particle  interaction,  the  sooner 
the  whole  bed  is  ignited.  The  particle-particle  interaction  law  used 
here  is  the  one  employed  by  Kuo  [22],  where  the  nominal  value 

K = 6.94  X lo"*  Ibf/in^  was  utilized.  Large  values  of  the  particle-par- 
ticle interaction  indicate  that  the  bed  would  resist  the  further  com- 
paction more  severely.  Correspondingly,  for  the  case  of  10  K,  apparently 
the  bed  resistance  to  compaction  is  so  high  that  before  the  complete  bed 
is  ignited  the  numerical  code  goes  unstable  due  to  resulting  excessively 
high  gradients.  The  flame  locus  does  not  show  any  interesting  flame 
accelerations  or  any  significant  change  in  the  total  bed  ignition  time. 

As  with  the  interphase  heat  and  momentum  correlations,  there  seems  to  be 
only  a limited  range  of  these  formulations  which  allow  (with  our  com- 
putational scheme)  for  well-behaved  solutions  which  lead  to  complete 
ignition  of  a given  length  of  granulated  propellant. 
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Figure  4.22  Flame  front  locus  as  a function  of  propellant  burning  rate  sensitive  to  the 
bulk  temperature  of  the  solid. 
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APPENDIX  A:  DERIVATION  OF  THE  TWO-PHASE 
CONSERVATION  EQUATIONS  WITH  SEPARATE-FLOW  CONCEPT 


The  following  is  a summary  of  the  key  steps  generally  used  to  derive 
the  six  field-balance  conservation  equations.  Some  of  the  details  were  pre- 
sented in  a graduate  course,  AAE  493,  taught  by  Professor  H.  Krier,  on  the 
subject  of  one -dimensional  two-phase  flow.  The  key  assumptions  are: 

1)  The  gas  phase  and  the  solid  phase  are  treated  separately. 

2)  The  motion  of  the  two  phases  is  coupled  through  interaction  terms 

3)  Each  phase  is  treated  as  a continuum. 

4)  Solid  particles  are  large  compared  to  molecules. 

Conservation  of  Mass 

(a)  Gas  Phase  Continuity 


Consider  a control  volume  with  gas  flowing  through  it. 


(p  u A ) 
g g g X 


(P  u A ) . 

g g g x+Ax 


increase  of  mass  inside  = mass  added  to  control  volume  in  time  t 


3(p  A Ax) 

= (puA)  -(puA)  . + r SAx  + r Ax 

g g g X ^"^g  g g"x+Ax  c w 


where 


p = density  of  the  gas 
g 

A = cross-sectional  area  of  gas  phase 
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Ug  = velocity  of  the  gas 

= gas  source  due  to  combustion  of  solid  particles,  = 
r = gas  source  due  to  combustion  of  walls 
S = total  cross-sectional  area  of  control  volume 
Note  that  F^  is  defined  as  mass  per  total  volume  per  unit  time.  Also  F^  is 
defined  as  mass  per  unit  length  per  unit  time. 

Now  divide  by  Ax  (which  is  not  a function  of  time)  and  let  Ax  0. 


3(p  A ) a(p  u A ) 

.g-g  = - -g-g-g-  + r S + F 

9t  9x  c^  w 


(la) 


Now  use  the  definition  of  porosity  to  further  simplify: 


A Ax  A 


A = (t)S 
g 


S-A  A 

1-^  = ^ 4 


A = (1-(|))S 

p 


Substituting  for  A 


g 


3(Pg0S) 

3t 


9(p  u 05) 
g.  g- 
3x 


F S + F 
c w 


(lb) 


Divide  by  S, 


3(0p  ) , 3(0pu  S)  F 

^ g-^  1 g _ p _w 

3t  S 3x  ■ c ■ S 


(Ic) 


Let  Pj  = 0Pg  and  expand  partial  derivatives. 


'fPl^  1 

at  s 


ac  Dp, 

p,u  + p,S  ^ + u S ^ 
1 g 3x  1 Dx  g ax 


= 'c^f 


(Id) 


Rearranging  the  term,  the  final  form  is  written  in  the  manner  most  useful 
for  implementation  in  a finite  difference  scheme,  that  is, 


I 

(2) 


(3) 


f 


[ 
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Note  that  F is  defined  as  mass  per  total  volume  per  unit  time.  Also  F _ 
p r r yip 

is  defined  as  mass  per  unit  length  per  unit  time.  Again  divide  by  Ax  and 
let  Ax  -►  0. 


3(p  A ) a(p  u A ) 

. r S . F 

ot  3x  p wp 


(3a) 


Substituting  for  A^  from  the  previous  page 


3[pu  (1-<^)S] 


at 


£_e: 


ax 


= F S + wp 
P 


(3b) 


Divide  by  S and  employ  the  relation  ‘ 


3I(1-'J>)P  ] 1 [(1-<I')P  u S]  F 

- P^  + i p p = . r + 

at  s ax  c s 


(3c) 


Let  = (l-(l))p^  and  expand  partial  derivatives 


3(P2)  1 

“5r”  ^ s 


ii 

*^2^p  ax  * ^2  ax 


ap. 

+ u S ^ 
p ax 


F 

= -r  * ^ 

c S 


(3d) 


Rearranging  terms  yields  the  final  form  of  the  equation 


(4) 
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Conservation  of  Momentum 
(a)  Gas  Phase  Balance 

Once  again  consider  the  control  volume  with  mass  flowing  through  it. 


Time  rate  of  increase 
of  momentum  inside 
control  volume. 


Thus 


The  sum  of  the  forces  acting 
on  the  control  volume.  These 
include : 

Momentum  flux  through  entrance 
and  exit. 

Interaction  forces. 

Pressure  gradient  force. 

Forces  associated  with  mass 
addition  or  loss. 


3(p  A u ) 

ct  rt  rt** 


3(p  A u^) 

~ rf  n or  * 


+ r u s - r u 

c p w wg 


A 

g 9x 


.Note  that  it  is  assumed  the  pressure  gradient  force  acts  only  on  the  area  A . 

S 

It  is  also  assumed  that  since  the  gas  is  formed  at  a velocity  equal  to  the 
particle  velocity,  there  is  a source  of  gas  momentum  per  unit  area  equal  to 


r u . 


c p 


Of  course  the  solid  phase  must  have  a sink  of  momentum  equal  to  that 


amount.  Also,  here. 
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u = X - component  of  the  velocity  of  the  gas  leaving  the 
control  volume. 

¥ ® viscous  interaction  force  between  the  gas  and  the  solid 

particles. 

To  simplify  Eq.  (5)  first  substitute  = (t)S: 


3(p„4>Su„  3(P„<t»Su2) 


+ rus-ru  -())S|^-FS 
c p w wg  ^ ax 


Divide  by  S and  let  (pp  = Pj^. 


a(PiU^)  . a(p,su|)  _ 

9t  S 9x  c p S 


9P  — 

- 35^-  F 


Now  expand  the  partial  derivatives: 


9u  9p,  , 9u  9p, 

Pi  " ^g  ^ " s [^Pi^^g  3^ " pi^g  n " % ir 

r u 

. ru  - JUS. . p |£  . F 

c p S 9x 

To  further  simplify  this  equation,  take  the  gas  continuity  equation  (2) 
and  multiply  it  by  u 


“g  sr-  -Vg  "g  Ik s 37*  “g'c  * s 


Subtract  Eq.  (2*)  from  Eq.  (5c)  to  get 


9u  9u 

D.  + O.u  -;r-^  = r (u  -U  ) - -^(u  +U  ) 

^1  9t  ^1  g cbc  o'-  p g-^  S '■  wg  g-^ 


F 

97  ■ ^ 


Rearranging  the  terms,  one  finally  obtains 
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(6) 


(b)  Solid  Phase  Momentum  Balance 

Remembering  the  assumption  that  the  particles  act  as  a continuum,  the 
solid  phase  momentum  equation  can  be  written  as  follows. 


5(p  A u ) 3(p  A u^)  3t  _ 

LP  P P^P  - r u S - A 3“^  + FS 

9t  3x  c p p 3x 


(7) 


Note  that  here  it  is  assumed  the  pressure  gradient  force  acts  only  on  the 

area  A and  particle  source  from  the  walls  has  been  neglected.  To  simplify 
P 

Eq.  (7),  first  substitute  A^  = (1-4>)S 


3rp  (l-(^)Su  ] a[p_(l-(|))Su2] 

js: 2_  E__ 2-  - r u s - (i-(|))s  3-^  + FS 


9t 


9x 


(7a) 


Divide  by  S and  let  (1-4))  Pp  = P2 


3(P2U„) 

3t  s ax 


r u 

c p 


(1-4)) 


(7b) 
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Now  expand  the  partial  derivatives: 

3u  , r 9u 

p + U ^ l2o  Sii 

^2  3t  p 3t  S 


3u  3p- 

2p-Su  T"  a — 

^2  p 3x  2 p 3x  p 9x 

3t 

= -ru  -(i-(|))  ■»  ■*  + F 
c p 3x 


C7c) 


To  further  simplify  Eq.  (7c),  multiply  the  solid  phase  continuity 


Eq.  (4)  by  Up. 


3p, 

^ 

■“p  3t 


u-  — = - P2“p  ^ • '^p  9^^ S H ■ ^c'^p 


(4*) 


Subtract  Eq.  (4*)  from  Eq.  (7c)  to  get 


3u 


3u 

^2%  3x 


2 9t  ^2“p  3x 


3t  _ 

- (l-*)  ^ * F 


(7d) 


Divide  by  P2  and  rearrange  terms: 


(8) 


Conservation  of  Energy 

(a)  Gas  Phase  Energy  Balance 

The  first  law  of  thermodynamics  for  a control  volume  provides  that 

PgVe  “ ^ 


- r s 

c 


gchem  ^ 
g 


u^ 

N 

P 

P A 

E * * 

U 

g g 

g 2 

"gj 

g 

(9) 
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Note:  1.  The  walls  are  neglected  as  a source  or  sink  of  gas, 

2.  E = /c  dT  = c^T  (calorically  perfect  gas) 

9 f 9t) 

3.  ~ ^ I conduction) 

4 gChem  _ combustion 

8 1^2 

5.  r^S  = kinetic  energy  added  because  gases  are  generated  with 

the  velocity  u^. 

6.  <I>  = interphase  energy  transfer  from  gas  to  particles  and  interphase 

dissipation  work 

= q + -^  {t  u (l-d) } 

^pg  3x  XX  g"- 

q = h(T  -T  ) (1-a) (S/V) 

Tg  ^ 8 P ^ ' ''p 


To  simplify  Eq.  (9)  , substitute  = (J>S  and  divide  by  S. 


9 ^ 

(J)p  q = Fu  + <I)  + 7^  (bp  E + 
g^e  p 3t  g 2j 


,13  sIe  . ^lu  -r 

s 3x  "^'^g  I g 2 PgJ  g eg  2 ' 


Let  Pj  = (J)Pg  and  expand  terms. 

3 f f 

p.q  *Fu  + 4>  + p,  -^E  + -^  + E + -^\  r— 

1 e p I3t|_g  2J  (g  2j3t 


u^  P 


u^  P 


1 pn  sf  E .p,u  E 

S 1 g 3x  g 2 p_  1 g g 2 p_  3x 


u^  P I 3u 


u:  PJ  9P, 


. PCE  ;:i.us  e 


g 2 PgJ  3x  V [ 8 2 ' PgJ  3x 


r 

■ eg  2 
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Divide  by  Pj  and  form  substantial  derivatives 


Fu  » ri  Si  ^ 

V ^ D_  g ^ g g 
pj  pj  Dti  g 2 J 8 Pg 


(E  + uV2)  p,  u P 9p, 

g g -IL  + _J. 

p,  Dt  p,p  9x 


r , 

C pChem  _g^ 

P,  \ 2 


u 

u^  P 1 

9S 

- — + 

u^  P 1 

+ -S. 

E + _g.  + -S. 

E + 

j 

s 

g 2 p 

9x 

In  order  to  further  simplify  Eq.  (9c) , multiply  the  gas  continuity  Eq.  (2) 

(E  - uV2) 
by  — ®— - — ® to  get 


E. 


p^  Dt 


un 

9u 

u^l 

u 

E + 

. 

E + — ^ 

_S. 

i g 2 J 

9x 

Ig  2j 

S 

9x 

CE„-^uV2D 


(Remember  it  is  now  assumed  that  = 0.) 
Now  add  Eq.  (2**)  to  Eq.  (9c) 

p^  pj  Dt  i g 2j  g - 


un 

9 

+ U — 

f \ 

P 

E + 

1 g 2j 

g 9x 

l^gJ 

r I u^  u^l  u P P 9u 

_c  chem  ^ _E  _ n _ + ..M.  S ^ ^ S.  S. 

Pi  I g 2 “^g  2 J SPg  9x  Pg  9x 


Now  multiply  the  gas  momentum  Eq.  (6)  by  u^  to  get; 


Du  u 41  r u u F 

u (U  - u ) - ^ 


g Dt 


Pj  3x 


g 


(6*) 


Add  Eq.  (6*)  to  Eq.  (9d)  and  rearrange  to  get: 


D(Ep  F(u„-u^) 
Dt  p 


£1  = . -i.  ^ q + 

P,  p, 


gChem  ^ _£ 

I g 2 , 


+ u u 
2 g p 


u P ~ P 8u 

Jg_g.  ^ 

Sp  3x  p 3x 

g g 


u P 
g g 9x 


l^gJ 


u P 3p 

PjPg  3x 


(9e) 


The  last  two  terras  can  be  combined  as  follows: 


^ f 1 u P 3p 

p--" 


g g 3x 


p2  3x 
g 


u P 3p,  u P (j)  3p  u P 
g g 1 ^ _ g g g _ g g 

p^pg  3x  (l)pgpg  3x  p^  3x 


Rewrite  Eq.  (9e) , substituting  for  q^  and  t. 
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(b)  Solid  Phase  Energy  Balance 

Expressed  in  a form  similar  to  the  gas  phase,  the  solid  phase  energy 


balance  is 


pA^  = -FuS-$S  + -^  p A (E  + uV2) 
P T P P P P P 


3 . r-  P g r c cChem  p 

3x  ^p  p pi  p 2 PpJ  c p 2 


Note  the  signs  on  the  terms  FUpS  and  <I>S  are  opposite  those  in  the  gas  phase 
energy  equation  because  the  overall  mixture  equation  (gas  + solid)  cannot  have 
such  source  terms.  For  this  derivation  the  following  form  for  the  particle 


enthalpy  has  been  assumed. 


u^  P 

H = E + ^ + 7^ 

P P 2 Pp 


That  is,  the  "continuum"  particles  act  like  a gas. 

To  simplify  Eq.  (11),  substitute  Ap  = (1-()))S  and  divide  by  S. 


Pp(l-<{»qp 


„ [ u^  P 

TT  P (l-4>)Su  E * 

S 3x  ^p^  P P 2 p 


p gChem  ^ 

c p 2 


(11a) 


Let  P2  = (l-<J>)Pp  and  expand  terms. 

^Generally,  would  be  a negative  number,  representing  the  latent  heat 

of  the  solid? 


r 2^ 

f .2'i 

u 

E > 

E ^ ^ 

[p  2 J 

IP  2 , 

ap, 

3t 


+ p-U  T— 

2 p 9x 


u''  P 

E + + -^ 

p 2 p 

^ P) 


S 


u-  P 1, 


E . ^ . -S.||S 
p 2 Ppj3x 


+ u 


- r 


P 

E 

P 2 p^i 


3x  " ^2 


u^  P ]3u 
E + ^ ^ 

[P  “ 


,chem 


. _E 

Ij 


Cllb) 


Rearrange  terms  to  get  substantial  derivatives. 


’2S 


D(E  +uV2) 

_ P P , 

f 2*' 

c + E. 

DP2 

Dt 

1 P 2 , 

Dt 

3 

+ P_U  -ir- 

[P  ' 

Li. 

^^2  i 

+ u - r 

f , u^ 

pChem  _£ 

^2  p 3x 

i^pJ 

p p 3x  c 

i P " 2 . 

.!2i. 

S 


r 


.2  p 


E . ^ 

P 2 Pp 


3S 

^ * P2 


u"  P )3u 


E * ^ > -S. 
P 2 Pp 


3x 


Cllc) 


To  further  simplify  this  equation  multiply  the  solid  phase  continuity 
Eq.  (4)  by  (Ep  + uV2). 


Dp,  9u  p u . 

H-  ■ - f2"p*''p/«5^  - 57  - ‘'c(Ep*u;/2) 


(4**3 


te  • • 


O 


Eq.  rilcl. 
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1 


^2^  9u 


3 

fP  1 

-S. 

+ u 

[p  1 

-S. 

[ap2  ] 

, ECh» 

3x 

.^p| 

p 

i^p] 

3x 

p p 

P-U  P ~c, 

+ 2 P & as 

Spp  8x 


Cl  Id) 


This  equation  can  be  further  simplified  by  multiplying  the  solid  phase 
momentum  Eq.  (8)  by  p2^p‘ 


Du  9t  _ 

PoU  ^ = - (l-<t))u^  + Fu^ 


'2“p  Dt  p ax 

Now  add  Eq.  (8*)  to  Eq.  (Hd). 


(8*) 


^2^ 


DE  p-P  au 

- 4)  + 0 ^ — 1 

^2  Dt  pp  ax 


u P 3p, 


p 8x 
P 


3t 
_i 
p 9x 


P P 


n PoU  P. 


sp^  ax 


(lie) 


3(E  +uV23 


f>2S  ■ - ’’“p  - % * ‘>2  — V- 


dp  3(E  +uV2) 

+ (E  ^u^/2)~  + p.u  

p p dt  2 p dx 


9u  (E  +uV2) 

+ P-(E^+uV2)  tt-E-  - p,u  — 7^ 


■'2'  p V ' 3x  ^2  p S 3x 


3p, 


+ u (E  +uV2)  - r - uV2] 

p'-  p p 3x  c*'  p p •' 


Once  again  simplify  by  adding  Eq.  (4**)  to  Eq.  (13a). 


DE  Du 

Pu  - $ + Po  TTr^  + P-,u  [E^  + E 


‘'2'^  " ■ ‘“p  ■ "p  ^ ^2  Dt~  ■ ^2“p  Dt  ‘c^'p 


Simplify  further  by  adding  Eq.  (8*)  to  Eq.  (13b). 


DE  3t  , 

L = - ‘J’  + Po  rvT^  - -5^  - r [E  + E ® 

^ p 2 Dt  p 3x  c*-  p p 
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A third  form  of  the  particle  phase  energy  equation  can  be  derived  from 
the  assumption  = E^.  With  this  assumption  a drop  in  the  particle  kinetic 
energy  cannot  change  the  particle  temperature.  The  solid  phase  energy 
balance  is 


P A Cl  = -FuS-$S+-5-rpAE] 
P P T P p 3t  c^^p  p pJ 


+ ^ [p  A u E ] - r S[e‘^’*^®"’  - uV2] 
P P P P c ^ p p'  ■> 


(15) 


Substitute  A^  = (l-(|))S,p2  = (l-<!))pp,  divide  by  S and  expand  terms. 


9E 


9p, 


9E 


* p ^ ^2Vp  9S  , ^P2 

" P2^p  ^ H " Vp  97" 


- ^ct^p*'""'  - ^p/2] 


(15a) 


Simplify  by  multiplying  the  particle  continuity  equation  by  E^. 

Dp-  9u  p u E 

E _i  = . o E —2.  - ^ . E r 

p Dt  ^2  p 9x  S 9x  pc 

Add  Eq.  (4***)  to  Eq.  (15a) 


DE 

p q^  = - Fu  - $ + p ^ - r + E - uV2] 

273  p p 2Dt  c p p p -* 


(15b) 


Divide  by  p^  and  rearrange  terms. 


8 


3E 

_B.  = 
9t 


3E 

- u 


p 9x 


Fu 

+ -£ 
P2 


^ _£.IE  + 

Po  P P 


u^/2] 


$ 

+ -£. 
Po 


(16) 


le  phase  momentum  equation  to  put  this 
era tor  form. 

the  particle  energy  equation  (Eqs.  12, 
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APPENDIX  B:  CONSTITUTIVE  RELATIONS 


As  stated  in  Chapter  Two  the  solution  of  the  conservation  equations 
requires  some  knowledge  of  the  physical  laws  that  govern  the  rate  and  inter- 
action processes  involved.  In  some  cases  use  of  these  laws  was  extrapolated 
to  areas  not  within  their  stated  domain.  The  applicability  of  some  of  these 
laws  is  discussed  by  Krier  and  Van  Tassel 1.^ 

A general  criterion  for  particle  ignition  assumes  ignition  to  have 

occurred  when  the  surface  temperature  of  a particle  reaches  a critical  value. 

22 

In  order  to  determine  the  surface  temperature,  Kuo  e^  al^.  solve  the  heat 
conduction  equation  for  the  particle  and  do  not  employ  a particle  energy  equa- 
tion. Although  this  approach  may  define  ignition  appropriately,  Krier^  argues 
strongly  that  any  separated  flow  model  must  satisfy  a mixture  energy  equation 
which  is  the  sum  of  separate  phase  energy  equations  (Eqs.  10,  12).  To  predict 
ignition  he  expresses  a solid  phase  energy  balance  (Eq.  11)  and  defines  a mean 
"bulk  temperature,"  as  being  an  average  temperature  of  the  solid.  Ignition 
occurs  when  this  "bulk  temperature"  reaches  a critical  value  called  which 

is  less  than  the  surface  critical  temperature.  This  latter  approach  was  taken 

here  and  a value  of  560°R  for  T,  was  used. 

ign 

The  convective  heat  transfer  between  the  hot  gases  and  the  solid  par- 

32 


tides  was  described  using  Denton's  formula 

K1 

l^p=  0.65  -S— 2 S E 


That  is 


K 


g 


p 0.  33 

r 


CB.l) 


There  is  some  question  whether  the  same  formula  should  be  used  for  inter -phase 
heat  transfer  before  and  after  ignition  of  the  particles.  It  seems  reasonable 
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to  believe  that  the  flow  field  around  burning  particles  differs  from  that 
around  inert  particles,  and  that  these  differences  act  to  change  (reduce) 
heat  transfer.  Because  of  this  effect  on  the  boundary  layer  by  the  genera- 
tion process,  the  above  heat  transfer  coefficient  was  reduced  after  ignition, 
usually  to  10%  of  its  actual  value. 

Two  laws  for  drag  between  the  particles  and  gases  were  used.  One  was 
Ergun's  correlation^^ 


-u„)  |u„-u„|  (1-(|)) 


150Cl-<}>) 


+ 1.75 


and  the  other  was  a correlation  devised  by  Kuo  and  Nydegger' 


_ y(u  -u  ) (1-4))^  Re 

^ ^ d^5 276.23  + 5.05 

V I 


CB.3) 


Here  the  Reynold's  number,  was  calculated  based  on  the  radius  of  a 

solid  particle. 


2r  p 4>  u - u 


Re  = 
P 


In  the  course  of  this  work  it  proved  desirable  to  reduce  both  of  these  drag 
laws  by  a constant  fraction.  Justification  for  this  reduction  is  based  on 
the  application  of  these  laws  to  regions  outside  their  stated  domain,  use  of 
three  dimensional  laws  to  a one -dimensional  problem,  and  uncertainty  in  the 
correct  gas  viscosity. 


For  the  burning  rate  law  the  simple  pressure  dependent  relation  was 


employed. 


r = b^  . b^CPg)- 


(B.5) 


Eq.  (B.2)  therefore  defines  f - = P/Cu  - u ) 

or'  o 
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The  production  rate  of  gases  from  solid  particles,  per  total  volume  of  the 
bed,  is  therefore  given  by 


I’g  = (1-^)  Vivj 


(l-<fr)Ppr 

“ d“  /6 

P 


(B.6) 


where 


is  the  surface  area  to  volume  ratio  for  a single  particle. 


The  expression  for  axial  normal  stress  between  particles  was  that  used 


by  Kuo  and  Summerfield,^^  that  is 


(l-(j))-MKr(l-(})J“'  - (1 


- 


if  ((>  - 4> 


if  4)  > <() 


CB.7) 


where  = .4  and  K = 6.94  x 10  Ibf/in^.  The  constant,  (})^,  is  a critical 
porosity  above  which  the  particles  do  not  touch. 

The  last  relation  of  concern  in  the  program  has  to  do  with  the  varia- 
tion of  viscosity  with  temperature.  Ideally,  a relation  using  both  tempera- 
ture and  pressure  should  be  used  since  with  the  extremely  high  pressure  en- 
countered it  is  natural  to  assume  that  viscosity  would  be  affected.  Unfor- 
tunately no  such  relation  could  be  found.  Two  expressions  were  employed  which 
vary  viscosity  as  a function  of  temperature  alone.  They  are 


and 


p = p. 


P = P 


1 . B5 


T 

_g. 

T 

o 


(T^/537) ^ ® 
Tfg+198) 


(B.8) 


where  p^  and  T^  are  viscosity  and  temperature  at  initial  conditions  and  p*  is 


not  a viscosity  but  a parameter  with  units  of  Ibm  - ®R/in-sec.  When 


82 


calculating  a heat  transfer  coefficient  between  gas  and  particles,  an 
average  temperature  K^g  * was  used  instead  of  just  in  the  second 

expression  to  calculate  an  appropriate  viscosity. 
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APPENDIX  C:  CALCULATION  OF  THE  SOUND  SPEED 
IN  THE  TWO -PHASE  MIXTURE 

In  order  to  calculate  the  sound  speed  of  the  mixture,  the  sound  speed 
in  each  phase  must  be  known. 

Convention:  1 - indicates  gas  phase 
2 - solid  phase 
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Now  the  sound  speed  of  the  mixture  can  be  calculated  using  the  following 
formula  from  Wallis  [8]. 


1 


ap^  + (l-a)Pj 


a 1-a 


= 9.46  X 10"^°  secVin^ 


C = 32,500  in/sec 
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